---
title: "`r config.population.name` stock assessment on data-limited CMSY model (`r year.start` - `r year.terminal`) in `r config.population.area`"
author: "Developer: Piatinskii M., Report build by: `r config.report.author`"
date: 'Report build date: `r Sys.time()`'
output: 
  html_document:
    toc: true
    toc_depth: 3
    toc_float: true
    theme: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = FALSE)
```

## 1. Model information
**CMSY** - The CMSY model developed by Froese et al. 2017 employs a stock reduction analysis using priors for r based on resilience, K based on maximum catch and the r priors, and start, intermediate, and final year saturation based on a set of simple rules. It also allows users to revise the default priors based on expert knowledge. The SRA employs a Schaefer biomass dynamics model and an algorithm for identifying feasible parameter combinations to estimate biomass, fishing mortality, and stock status (i.e., B/BMSY, F/FMSY) time series and biological/management quantities (i.e., r, K, MSY, BMSY, FMSY).

To perform CMSY model only catch information required. 

Reference: [Froese R, Demirel N, Coro G, Kleisner KM, Winker H (2017) Estimating fisheries reference points from catch and resilience. Fish and Fisheries 18(3): 506-526.](http://onlinelibrary.wiley.com/doi/10.1111/faf.12190/abstract)


## 2. Input data
This section present the input data for stock assessment procedure. 

Required input data:

  - C - catch information over years. Skips not allowed.

```{r}
knitr::kable(data, caption = "Table 2.1. Input data", align = "l")
```

## 3. Model tuning
Model resilience tuning: **`r config.population.resilience`**

CMSY model performed using next cmd: 

```{r eval = FALSE, echo = TRUE}
cmsy <- cmsy2(year=data$year, catch=data$catch, resilience = config.population.resilience)
```

## 4. Results
There summary CMSY modelling results shown. Feel free to use it.

### 4.1 Estimates - B, F

```{r echo = FALSE}
res <- cmsy$ref_ts %>%
  select(year, b, b_lo, b_hi, f, f_lo, f_hi) %>%
  mutate_at(c(3,4), ~round(., 0)) %>%
  mutate_at(c(6,7), ~round(., 3)) %>%
  unite("b.ci95", b_lo:b_hi, sep = " - ", remove = TRUE) %>%
  unite("f.ci95", f_lo:f_hi, sep=" - ", remove = TRUE)

knitr::kable(res, caption = "Table 4.1.1. Biomass and fishing mortality estimates", align="l")
```

Column description:

  - b - biomass estimation
  - b.ci95 - biomass confidence interval at p = 0.95 level
  - f - fishing mortality estimation
  - f.ci95 - fishing mortality conf.interval at p = 0.95
  
```{r echo = FALSE, fig.cap="Fig. 4.1.1. B estimates with 95% confidence interval"}
plot(cmsy$ref_ts$year, cmsy$ref_ts$b, xlab="Year", ylab="B", type="b", ylim = c(0,max(cmsy$ref_ts$b_hi)))
lines(cmsy$ref_ts$year, cmsy$ref_ts$b_lo, col="red", lty=2)
lines(cmsy$ref_ts$year, cmsy$ref_ts$b_hi, col="red", lty=2)
```

```{r echo = FALSE, fig.cap="Fig. 4.1.2. F estimates with 95% confidence interval"}
plot(cmsy$ref_ts$year, cmsy$ref_ts$f, xlab="Year", ylab="F", type="b", ylim = c(0,max(cmsy$ref_ts$f_hi)))
lines(cmsy$ref_ts$year, cmsy$ref_ts$f_lo, col="red", lty=2)
lines(cmsy$ref_ts$year, cmsy$ref_ts$f_hi, col="red", lty=2)
```
  
### 4.2 Reference points
MSY reference point stategy done. Mean reference points shown in table 4.2.1.

```{r echo=FALSE}
  cmsy$ref_pts[-2:-1,] %>%
  mutate_if(is.numeric, ~round(., 3)) %>%
  mutate_if(is.numeric, ~format(., scientific = FALSE)) %>%
  knitr::kable(., caption="Table 4.2.1. MSY reference point estimates", align="l", row.names=FALSE, digits=4)
```

Column description:

  - param - reference point name column
  - est - estimated value
  - lo - lower confidence level
  - hi - upper (higher) confidence level

Remember, Bmsy calculated once (one time) in whole time seria. Fmsy point calculated for every one year but not presented here.

### 4.3 Estimates VS reference points

```{r echo = FALSE, fig.cap="Fig. 4.3.1. B/Bmsy proportion in time vector"}
plot(cmsy$ref_ts$year, cmsy$ref_ts$bbmsy, xlab="Year", ylab="B/Bmsy", type="b", ylim = c(0, 1.5))
lines(cmsy$ref_ts$year, cmsy$ref_ts$bbmsy_lo, col="red", lty=2)
lines(cmsy$ref_ts$year, cmsy$ref_ts$bbmsy_hi, col="red", lty=2)
abline(h=1, col="blue", lty=3)
```

```{r echo = FALSE, fig.cap="Fig. 4.3.2. F/Fmsy proportion in time vector"}
plot(cmsy$ref_ts$year, cmsy$ref_ts$ffmsy, xlab="Year", ylab="F/Fmsy", type="b", ylim = c(0, 3))
lines(cmsy$ref_ts$year, cmsy$ref_ts$ffmsy_lo, col="red", lty=2)
lines(cmsy$ref_ts$year, cmsy$ref_ts$ffmsy_hi, col="red", lty=2)
abline(h=1, col="blue", lty=3)
```

```{r echo = FALSE, fig.cap="Fig. 4.3.3. Catch, biomass, MSY and Bmsy in time vector"}
plot(cmsy$ref_ts$year, cmsy$ref_ts$b, xlab="Year", ylab="Biomass/Catch", type="b", col = "blue", ylim = c(0, b.max))
abline(h = cmsy$ref_pts$est[3], lty = "dashed", col = "red")
lines(cmsy$ref_ts$year, cmsy$ref_ts$catch, col="red", lty=2, type = "b")
abline(h = cmsy$ref_pts$est[5], lty = "dashed", col = "blue")
legend("topleft", legend = c("B & Bmsy", "C & MSY"), col = c("blue", "red"), lty=c(1,1))
```

### 4.4 Model r/K parametrization
```{r echo=FALSE}
  cmsy$ref_pts[1:2,] %>%
  mutate_if(is.numeric, ~round(., 3)) %>%
  mutate_if(is.numeric, ~format(., scientific = FALSE)) %>%
  knitr::kable(., caption="Table 4.4.1. Schaefer model r/K optimum found", align="l", row.names=FALSE, digits=4)
```

Column description:

  - param - param name column
  - est - estimated value
  - lo - lower confidence level
  - hi - upper (higher) confidence level

## 5. Diagnostics
Some avaiable model diagnostics are summarized here.

### 5.1. Retrospective analysis
There is summary retrospective figure presented. Retrospective procedure are syntetic continious reducing time seria by 1 year.

```{r echo=FALSE, fig.cap="Fig. 5.1. Retrospective B analysis"}
plot(cmsy$ref_ts$year, cmsy$ref_ts$b, type="l", ylim = c(0, b.max), lwd=3, xlab="Year", ylab = "B, biomass")
for (i in 1:config.retro.years) {
  lines(retro.data[[i]]$year, retro.data[[i]]$b, col=(i+1))
}

```


```{r echo=FALSE, fig.cap="Fig. 5.2. Retrospective F analysis"}
plot(cmsy$ref_ts$year, cmsy$ref_ts$f, type="l", lwd=3, ylim=c(0, 2), xlab="Year", ylab="F")
for (i in 1:config.retro.years) {
  lines(retro.data[[i]]$year, retro.data[[i]]$f, col=(i+1))
}
for (i in seq(0, 2, by=0.5)) {
  abline(h=i, lty=2)
}
```

### 5.2. Mohn-rho test
The basic ICES (2018) procedure to determine model stability over years - retrospective Mohn rho test. Mohn's rho test calculate relative bias for latest 5 years (default ICES procedure - 5 years) retrospective variations on scale -1 ... +1. Low negative values of *rho* leads to underestimate factor, high positive values - to overestimation (remark: values lowest -0.4 or higher +0.4 shows high variations and low model stability). So procedure is:

$$
\begin{aligned}
relbias = (retro - base) / base \\
rho = mean(relbias)
\end{aligned}
$$

For SSB and F math approach will be

$$
\begin{aligned}
\rho_{SSB} = \frac{1}{n} \sum_{i=-5}^{n=0} \frac{(B_{i} - \overline{B})}{\overline{B}} \\
\rho_{Fbar} = \frac{1}{n} \sum_{i=-5}^{n=0} \frac{(F_{i} - \overline{F})}{\overline{F}} \\
\end{aligned}
$$

where i - year steps from terminal year (terminal-5, terminal-4, ... terminal).

```{r echo=FALSE}
retro.rho.ssbdata %>%
  rownames_to_column(., var="year") %>%
  mutate_if(is.numeric, ~round(.,1)) %>%
  knitr::kable(., caption = "Table 5.2.1. B retrospective values", row.names = TRUE)
```

```{r echo=FALSE}
retro.rho.fdata %>%
  rownames_to_column(., var="year") %>%
  mutate_if(is.numeric, ~round(.,3)) %>%
  knitr::kable(., caption = "Table 5.2.2. F retrospective values", row.names = TRUE)
```

**Final** Rho estimates: 

$$\rho_{SSB} = `r round(retro.rho.ssb,3)`$$
$$\rho_{Fbar} = `r round(retro.rho.f,3)`$$

## 6. Summary
There is summary plot from datalimited2 package displayed. This figures can be used directly to fishing regulation management.

```{r echo=FALSE}
plot_dlm(cmsy)
```

## Appendix list

### Appendix 1. Model reference points

```{r echo=FALSE}
print(cmsy$ref_pts)
```

### Appendix 2. Model estimates
```{r echo=FALSE}
print(cmsy$ref_ts)
```